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ABSTRACT 


This thesis analyses data of and Builds a simulation 
model for the track “et an underwater vehicle as perceived 
by a test range of three dimensional short baseline sonar 
arrays. In this way many random replications of track 
become avallable quickly and inexpensively. These 
Simulations support a larger project whose object is to 
monitor the performance of the test range and provide clues 
for troubleshooting problems. In particular, jyoint values 
of sensor array estimated displacement and reorientation 
corrections are generated and their error distribution is 


quantified. 
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I. INTRODUCTION 


The Naval Undersea Weapons Engineering Station (NUWES) 
operates nine test ranges that serve a variety of purposes 
for testing and validating the Navy’s current and future 
weapons systems. This thesis is in support of a larger 
project whose goal is to monitor the performance of the 
short baseline ranges (ie. those ranges in which each array 
produces a three dimensional track of moving vehicles). 
More specifically, it deals with the development of a 
computer Simulation which generates realistic random 
replications of an underwater vehicle’s track (ie. vehicle 
position perceived by a sensor array on the range at a 
number of equally spaced time points). Such simulations 
serve to provide information about the inherent variability 
in the output of the tracking range, especially for 
assessing the calibration error. Its use can provide an 
improved understanding of the processes involved and 4a 
reduction in the frequency of range shutdown for the purpose 
of resurveying the remote sensor locations. 

When the simulated replicate tracks are passed to the 
program KEKEYMAIN (a project FORTRAN program that estimates 
displacement and orientation corrections to the range remote 
sensor arrays), the output is used to generate bDivariate 


scatter plots of these corrections. Although extensive work 


of this type was not possible in the current thesis, the 
limited results indicate a surprising amount of variability 
and suggest that the nature and extent of variability can 
change noticeably with relatively minor changes in the 
localized conditions. Considerable testing using this 
simulation tool is clearly indicated. 

The research reported here involves three general 
activities: 


(a) Analysis of real underwater track data to learn their 
important behavioral characteristics. 


(6b) Simulation model formulation and construction. 

(c) Development and programming of algorithms to merge 
with the existing project programs to produce. any 
number of random replications with correction 
estimates for each simulated track. 

The organization of the thesis is as follows: 
Section Il contains explicit background material and 
provides a framework for the research and shows how it 
relates to the larger project. Section III explains the 
data analysis procedures and results. The simulation model 
is developed in Section IV and the interfacing with the 
Overall project programs is described in Section V. Results 
and conclusions are detailed in Section VI, with areas for 
future work listed in Section VII. 

A number of appendices are included to provide the 
detailed support too extensive to be included in the main 


body of the paper. They are referenced in the appropriate 


sections. Additionally, to speed the completion and 


availability of the KEYMAIN program, the author wrote two 
subroutines (CONECT and REDUCE) to be included in that 
package. Appendices E through G contain the development of 


these subroutines and the FORTRAN code listings. 
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Il. BACKGROUND 


Consider a three dimensional underwater tracking array 
composed of four hydrophones arranged tao define a loacal 
Cartesian coordinate system (see Figure 1). Such an array 
is capable of tracking a target downrange, crossrange and in 
depth in its local coordinate system; it is termed a short 
baseline array because of the short distance of the xX. /Y, 
and Z hydrophones from the "corner" hydrophone (30 feet) 
that is used to gather the three dimensional tracking 
information. 

The arrays take "fixes" of target position. A fix 15 
defined by a bearing (azimuth and elevation) and distance 
from the array to the target. Tracked targets on the range 
are fitted with an acoustical device, called a "pinger", 
that emits pulses of sound, or "pings", at a specific 
frequency at precisely timed, regular intervals. The 
elapsed time from the generation of the ping at the target 
to the reception of the ping at the array establishes a 
distance from the target to the array. Because of the 
distance that separates the individual hydrophones of the 
array, the ping arrives at each hydrophone at a slightly 
different time. This time difference can be resolved into a 
bearing from the array to target. (Actually, these fixes 


are “apparent” fixes. They are used to initialize a sound 
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ray tracing algorithm that leads to the "real" fix.) - 
series of such fixes establishes a target track in that 
particular array’s own local Cartesian coordinate system. 
If the precise position of the array is known, its local 
track information can be translated to one common range 
coordinate system. 

Position of an individual array in the range Cartesian 
coordinate system is described not only in terms of its 
downrange, crossrange and depth (X,Y and 2Z) # £coordinates 
(termed location), but also with respect to X-tilt, Y-tilt 
and Z-rotation (commonly called roll, pitch and yaw) angles 
from the range coordinate system axes (termed orientation). 
Both sets of measures are needed to translate accurately 
from local to range coordinates. 

Typically, a range 1s composed af many individual 
arrays. Nanoose Range, for example, has 24 arrays, while 
Dabob Bay has 7. The arrays are arranged in such a way that 
the array coverages overlap one another to provide 
continuous tracking on a target vehicle. Such overlapping 
areas are called crossover regions (see Figure 2), and 
produce two sets of track on the same target for the same 
time period. 

Ideally, the corresponding points, or fixes described 
earlier, from each array inaéecroassover region should 
translate to identical points on the range coordinate 


system. In practice, this seldom happens. 
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Three major sources of this variability in crossover 
track data are: 


(a) Slippage of the sensor arrays from their assumed 
positions in the range coordinate system. 


(b) Time synchronization and instrumentation problems. 


(c) Inhomogeneities and temporal variability in the 
water column. 


Although the reasons for slippage of the arrays are 
speculative, the fact that slippage occurs is evidenced by 
the change in sensor locations after a range resurvey is 
performed. The question of discriminating timing errors, 
item (b), from slippage errors is on a ‘future agenda. 
Investigation dane by Main (CRef. 1] provides us with 
methodology for treating the random components of these 
errors. Item (c) falls into a broader category which may 
require extensive investigation. It may also provoke a 
reassessment of the assumed uniform horizontal quality of 
the range water column. It seems wise to account for 
Slippage and instrumentation problems first. Our work 1s 
confined to (a). 

The present research is in direct support of Protessor 
Robert R. Read of the Naval Fostgraduate School who has been 
working on the slippage question. He has produced a FORTRAN 
program (KEYMAIN) which takes the crossover data and 
estimates array position corrections using a non 1inear 
least squares algorithm. The surfaces that enter into this 


least squares optimization function are rather flat and the 
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values do not change much as the location corrections are 
varied. Replicated data is needed to quantify the extent to 
which the estimated corrections are within the range of 
natural variability. The simulation model in this’ thesis 
provides the tool by which this’) variability can be 
quantified. With the simulation, the needed replicated data 
can be quickly and easily generated, and scatterplots of 
displacement (change in location) versus rotation (change in 
orientation) can be made to investigate the inherent 


variability in the correction estimates. 
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IItT. DATA ANALYSIS 


A. GENERATING RESIDUALS FROM REAL TRACK DATA 

The real track data used for this study was taken from 
Nanoose range in September, 1982 (Figure 3). It was 
supplied as an N x 7 matrix, N representing the number. of 
data points in the crossover data set and the columns 
representing the following: 

Column 1 Point count. Every pulse emitted by the 
Pinger is assigned a sequential number 
beginning at the start of the tracking run. 
Therefore, this column’s value increases by 
One each row unless a data point is missing, 
which would be indicated by a sequential 
omission in this first column. 

Columns 2-4 The X,Y, and Z coordinates of the target 
for the column one point count as determined 
by the first sensor in the crossover pair. 

Columns 3-7 The X,Y, and Z coordinates of the target for 
the column one point count as determined by 
the second sensor in the crossover pair. 

There were six such data sets used. These particular 
sets were chosen because each formed a straight line in the 
range space and so the best straight line path of the target 
could be estimated from the data. (In contrast, estimating 
the best curved path from curved track data would have been 
far more challenging, but with little or no gain in the 
examination of residuals.) 


The idea was to take the track data, fit the best 


straight line possible through the data, then examine the 
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residuals formed by subtracting the fitted straight line 
Point sequences from the original data. 

The commented FORTRAN program written to generate the 
residuals is included as Appendix A, and its mechanics are 
discussed below. 

Since the data was Given in a esingle matrix and 
residuals were desired for each sensor, the data was first 
separated into two tracks, each an N x 3 matrix; call them 
TRI and TR2. Next, a 3 x» 3 covariance matrix for each track 


was computed by the formula 


N 
COV = €)) (TRi)-TR)’ (TR(i)-TR)I] / N-1 


i=1 
where 
TR(i) = 3 component row vector of the X,Y,Z coordinates of 
the ‘eal data point 
AR = = component row vector of the average of the N 
(X,Y,2Z) coordinates for that track data 
N = number of data points 


The best straight line estimate of the target path 15 
the straight line such that the sum of the distances of each 
data point to the line 1S a minimum. This implies that the 
distance from each point to the line 1S a minimum - that 15s, 
each point should be projected onto the line orthogonally. 
The method employed to accomplish this orthogonal regression 
was that of principle components. Principle components 


requires that the eigenvector associated with the largest 


cy 


eigenvalue of the covariance matrix be identified. It was 
computationally convenient to make a 3 x 3 matrix of the 
eigenvectors, with the first column being the eigenvector 
associated with the largest eigenvalue, the second column, 
the elgenvector associated with the second largest 
eigenvalue and third column, the eigenvector associated with 
the smallest eigenvalue. This done, the following matrix 
multiplicatian 


PROJECTION = (TR —- AVE) x EIGENVECTORS 


where 
TR = N x 3 track data matrix 
AVE = N x 3S matrix made of N identical rows of 
the X,Y, and Z averages of the track data 
EIGENVECTORS = 3 X 3 eigenvector matrix described above 


yields the N x 3 matrix PROJECTION, whose entries are the 
orthogonal projections oaf the track data points onto the 
axes Of a new 3 dimensional Cartesian coordinate system 
whose origin is located at TR (the (X,Y,Z) averages for the 
data set) and whose axes are rotated such that the new xX 
axis is the best straight line that describes target path. 
For example, the first row of PROJECTION is a three 
component row vector; call the components (pl1l,p2,p3). Then 
the point (p1,0,0) is the projection of the first track data 
point anto the new X axis (known to be the best straight 


line that describes the target path), the point (0,p2,0) is 


the projection of the first track data point onto the new YY 
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axis, and the point (0,0,p3) is the projection of the first 
track data point onto the new Z axis. Since we are only 
interested in the projections onto the new X axis, we can 
replace the second and third columns of PROJECTION with 
zeros and have the coordinates of the best straight line 
path of the target. These coordinates are in the PROJECTION 
coordinate system, and must be translated back into the 
range coordinate system. Using the matrix PROJECTION (with 
last two columns = Q) and ~operforming the following 
multiplication 

OLD = EIGENVECTORS x PROJECTION’ 
yields the 3 x N matrix of the track data points in the 
range coordinate system. Note that OLD must be transposed 
to get the N x 3 format desired. Since the mean had been 
subtracted off before computing the PROJECTION matrix, to 
get the data points back to the proper position requires 
that the means be added back in. The matrix addition 

OLD’ + AVE 

yields the N x 3 matrix of the orthogonal projections of the 
data onto the straight line path of the target in the 
original range coordinate system. 

Actually, after having computed the PROJECTION matrix, 
there is yet another step to take before translating back 
into the range coordinate system. Based on an assumption of 
constant speed for the target, the track data points should 


be equally spaced since the pulses are emitted at regular, 


precisely timed intervals. This could be accomplished by 
using an interval equal to the total distance divided by the 
number of data points, but if data points are missing in the 
track aaee set (and each set of track data was missing 
several points), then this method introduces an error. The 
method chosen should show that, if a single data point is 
missing, the gap between the consecutive points on the data 
track should be twice as long as if none were missing. The 
first column (¢ the point counts) of the N x 7 track data 
matrix contains the information on missing points. 
Performing a simple least Squares regression of the point 
counts onto the first column of the FROJECTION 


matrixtranslates this information into the best straight 


line point sequence for target motion. Now when the 
projections onto the first principle component are 
translated back onto the range coordinate system, the 


result is truly the best set of track points attainable 
using all the information contained in the track data set. 
To obtain the residuals, the straight line estimated 
track points are subtracted, point for point, from the 
Original track data. It was important to discover the 
distribution of the residuals for each track segment because 
to simulate the track segments later, residuals were 
simulated and then added to the straight line. Knowing the 
character of the residuals allowed the generation of very 


realistic track data for the simulation model. 


B. ANALYZING TRACK DATA RESIDUALS 

Fach of the six track segments produced two sets of 
residuals - one set for each of the two sensors of the 
crossover pair. Each set was further divided into X, Y, and 
Z components. There were, therefore, 36 sets of residuals to 
be examined for distributional characteristics. Graphical 
analysis was performed on each set of residuals (histograms, 
cumulative density plots and Q@Q@ plots) and, when, from that 
analysis, a distribution for the residuals could be 
determined, formal statistical tests were performed to 
verify that the best distribution was chosen. The analysis 
Showed the residuals to be normally distributed. The best 
and worst case graphical and analytical results are included 
in Appendix &. Note that there were some very good fits to 
a normal density (pp. 66-71) with statistical significances 
for the Chi Squared and Kolmogorov- Smirnov goodness of fit 


tests well above avery conservative .35 1n all But one 


case. Even in the data that most poorly resembled a normal 
density (pp. 72-77), the Kolmogorov —- Smirnov goodness of 
fit significance level never fell below .25. From this 


analysis it was concluded that a normal density for the 
residuals was accurate and appropriate for the simulation 


model. 


C. TIME SERIES ANALYSIS 
From a simulation point of view, it 1s very convenient 
to assume independence of the residuals between successive 


2S 


points in the original track. If the assumption is not 
made, then one must resort either to using a point count 
interval of some integer value greater than one for making 
the simulated track from the sritesiaaie or Build = an 
autoregressive model. 

Concern for this correlation of the data over time led 
to a time series analysis of the residuals. This was done 
in two ways for each track from a single sensor. Recall 
that for each data point, residuals were generated in the X, 
Y, and Z directions of the range coordinate system. Each of 
these components was subjected to a time series analysis to 
determine whether there was a dependence in the errors in 
any single direction over time. Then, the distance of the 
data point from the straight line target path, given by 

SGRTC (RES *) + (RES. *) + (RES) J 
was examined to determine if the magnitude of the error of 
one point was correlated with the magnitude of the error of 
the next point. The time series analysis examines’ the 
dependence from point to point, on every other point, on 
every third point, and so on, so that the interval at which 
one can assume independence (indicated by an autocovariance 
Value of zero) can be determined. Appendix C contains’) the 
autocorrelation graphs for the & data sets; first the error 
magnitudes are examined for both sensors of a data set, 
followed by the individual error analysis in the X, Y, and Z 


directions. 
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There was no clear cut interval for which independence 


seemed to emerge in the data sets. There is instead a 
random pattern of insignificant dependence fram the 
beginning of the autocorrelation analysis, leading to the 


conclusion that point to point independence was appropriate 
for the simulation model. Attention 1s directed to the 
noise in the autocorrelation graphs with no distinct 
patterns emerging, and correlation magnitudes smaller than 
0.2 in most cases. It was decided, therefore, that any 
important time correlation was not so large as to cause 


excessive or even noticeable error in the simulated tracks. 


D. STUDY OF RESIDUALS 

The study of the residuals yielded important interesting 
information, Summarized in Table 1. The first column of 
Table 11s the standard deviation of the residuals from the 
left sensor in the downrange (X), crossrange (Y) and depth 
(Z) directions, three values for each data set. Column 8 
Gives the corresponding information for the right sensor oaf 
the crossover data pair. Columns 2, 3 and 4 for each data 
set are the columns of a3 x 3 correlation matrix for the 
residuals of the left sensor, columns 95, 6&6  @nd 7, 
corresponding information for the right sensor. 

First, note the difference between the left and right 
sensor standard deviations of residuals in all three 
directions. Although there are some very good comparisons 
(data set D2A1l, crossrange, and data set DSA, downrange) 
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TABLE 1 


STANDARD DEVIATIONS OF RESIDUALS 
AND 


CORRELATION COEFFICIENTS OF RESIDUALS 


CORRELATION 
LEFT 

DATA SET 

1.000 0.916 -—-0.899 

0.9186 Pa DOO S—On 7 2a 

-0.899 -0.728 1.000 
DATA SET 

1.000 oO; S355 "—O. 7a7 

0.859 1.9000 -0.7359 

-0.787 =-—-0.759 1.000 
DATA SET 

1.000 ©o.788 -0.888 

0.788 1.900 --0.,434 

-0.888 -—-0.434 1.000 
DATA SET 

1.9000 O.741 —-0.904 

Oey 41 1.000 -0.687 

-—0O.9704 -—-0.687 1.900 
DATA SET 

1.000 -0.930 -0.,.429 

—0. 920 1.000 re) Se 
—-0.429 a ye 1.000 
DATA SET 

1.000 -0.408 -0.450 

-0. 4093 1.000 -0. 454 
-0.4530 —-0,454 1.000 


1.000 


1.040 
Orel 
1.9000 


0.902 
-0.728 
1.000 


O. 390 
—-0.146 
1.900 


MO. 16e8 
-0.148 
12008 
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CORRELATION 
RIGHT 
D2A1 , N = 6&7 
1.900 -0.843 
-O.8635 1.000 
0.787 -9.455 
D2A2 , N = 70 
1.9000 —-0.260 
—-0. 3560 17 Ooo 
0.040 —-0.547 
D2B , N = S35 
1.000 -0.949 
—0.949 1.9000 
0.902 --0.728 
D4 s N= 93 
SWIG, Oa es 
10) 5 aye0 Pe Oe 
O.580 —-0.146 
DSA = Se 
130090 - Boe 59 oS 
104 seo 1.000 
O.168 -—-0.148 
DSB »5 N = 87 
1.608 0.445 
0.445 eee, 
O29 =O. 27S 
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there are also some wide disparities (data set D2B5, 
crossrange and data set DSB, crossrange). Using the 
variance ratio test described in Larson (Ref. 2: pp. 449], a 
statistical test was performed to determine whether the 
variances (the recorded standard deviations, squared) of the 
left and right sensor arrays for downrange, crossrange and 
depth were the same for a given data set. The results are 
given in Table 2. Using 90.05 as the basis for rejection of 
the null hypothesis (Ho: the two variances are the same), 8 
of the 18, or 44%, of the individual tests fail. With a 
confidence level of .95, one cate expect roughly i failure 
in the 18 tests. Eight failures is strong evidence that the 
standard deviation figures do not match very well. 

Recalling the test data from Figure 3 (p. 18), at 48s 
seen that there are three sensor arrays that contributed the 
data: arrays 4, 3S and 6. In data sets D2A1, D2A2 and D2B, 
array 4 18s the left array and array 3 is the right sensor 
array. For data sets D4, DSA and DSB, sensor array S is the 
left array and array 6 1s the right array. Ge waigne expect 
that the standard deviations in any single direction would 
be the same for a given sensor. This turned out not to be 
the case. Using Bartlett’s test for the equality of several 
variances (Ref. 3: pp. 225-2271, the hypothesis that the 
variances for a given sensor in a given direction are equal 
was tested. The results are given in Table 3. The first 


two sections of the table compare the 3 standard deviation 
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TABLE 2 TEST OF THE HYPOTHESIS THAT THE 
STANDARD DEVIATIONS WITHIN A DATA SET 


ARE EQUAL 
VARIANCE DEGREES LEVEL ADE 
RATIO OF FREEDOM SIGNIFICANCE 

DATA SET D2A1 

DOWNRANGE 1.363 66 -211 

CROSSRANGE - 788 _ 66 F935 

DEPTH 1.363 66 222? 
DATA SET D2A2 

DOWNRANGE 1.627 69 29045 

CROSSRANGE 1.244 69 2 564 

DEPTH 1.858 69 -O11 
DATA SET D2B 

DOWNRANGE 347 v2 -O31 

CROSSRANGE 2So7 w2 -90003 

DEPTH -378 v2 067 
DATA SET D4 

DOWNRANGE . 782 92 . 884 

CROSSRANGE 628 72 a Oa 4 

DEPTH 2692 92 979 
DATA SET DSA 

DOWNRANGE 2793 81 2939 

CROSSRANGE . 749 81 2195 

DEPTH . 908 81 -003 
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TABLE 2 


DATA SET DSB 


DOWNRANGE 1.527 
CROSSRANGE 3.300 
DEPTH - 638 


(continued) 
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Table 3 TEST OF THE HYPOTHESIS THAT THE 
STANDARD DEVIATIONS OF RESIDUALS OF 
A PARTICULAR SENSOR ARE EQUAL 


Chi Squared Random Variables 
2999 QUANTILE 


DEGREES OF FREEDOM 299 QUANTILE 


2 9.21 13.82 
= 15.09 20.52 
DOWN CROSS DEGREES 
RANGE RANGE DEPTH OF FREEDOM 
SENSOR 4 (LEFT) 39.57 2.02 1.41 2 
SENSOR 3S (RIGHT) 14.75 16.25 13.96 2 
SENSOR S (LEFT) 131.14 101.33 91.40 2 
SENSOR 6 (RIGHT) 148.72 175.18 B4.31 2 
SENSOR S&S (LEFT AND 153.04 150.87 105.951 ~ 
RIGHT) 

LEFT SENSORS 171.45 162.29 94.48 r=) 
RIGHT SENSORS 168.37 2502.97 107.15 2 


values for each sensor ina single direction. The test 
statistic is distributed as a Chi Squared random variable 
with 2 degrees of freedom. The .99 and .999 quantiles for 
such a random variable are 9.21 and 13.82. In only two 
cases of the 12 trials would the 3 standard deviations be 
considered statistically the same. 

“The third section of Table 3 is included because sensor 
array S was used as both a right and left array. Therefore, 


there were actually six values for downrange, crossrange and 


depth for that particular array. Using the same test as 
before, the question of whether the six values were 
statistically the same was’ investigated. The tabulated 


statistic for the 3 cases is distributed as a Chi Squared 
random variable with 5S degrees of freedom. The .99 and .999 
quantiles for such a ramdom variable are 15.09 and 20.52. 
In every case, the p-value of the test is essentially equal 
to 0.0. The results are highly significant. 

Finally, the bottom of the table investigates whether 
downrange, crossrange and depth residual standard deviations 
are the same for left and right sensors in general. 
Employing the same test for the 6 values produced the 
tabulated statistics. Since the test statistic 1S again Chi 
Squared with 5S degrees of freedom, the level of significance 
in every case is essentially 0O.O. 

A third question of interest from Table 1! is whether 


the correlation coefficients are the same for the left and 


—_ 


Sl 


right sensors ina data set. This was investigated in the 
following manner. 

First, the normalizing, inverse hyperbolic tangent 
transformation (Ref. 3: p. 3465] was applied to each of the 
off diagonal correlation coefficeients in the correlation 
matrices. (Note that the matrices are symmetric, so there 
are only three values of concern for each matrix.) This 


transformation makes each of the values normal with mean 





1 1 + rho 
- )\7%Fe—oeor———_—> 
2 1 - rho 
and variance 
1 
N—- 3 


where N is the number of data points contributing to the 
correlation. Under the assumption that the two independent 
sample correlation coefficients come *¢from the same 
population, their difference is distributed normally with 
mean 9 and variance 


es 





Ne > SS 
We can therefore go to the standard normal tables to obtain 


the significance of the statistic 


SORT C2/ (N-S) J 
which tests the hypothesis that the two correlation 
coefficients are equal. These values are tabulated in Table 


oz 


TABLE 4 TEST OF THE HYPOTHESIS THAT THE 
CORRELATION COEFFICIENTS WITHIN A DATA SET 


ARE EQUAL 
DATA CORRELATION STANDARDIZED LEVEL OF 
SET COEFFICIENT Zi - 22 NORMAL SIGNIFICANCE 

D2A1 X,Y 2.867 16.227 ond 
D214 ey 2.534 -14.333 0.0 
D2A1 Y,Z -.434 -2.453 0.014 
D2A2 X,Y 1.668 9.654 0.0 
D2A2 Kaa 1.103 -6. 384 0.0 
D2A2 Y,Z —. 381 -2.203 0.028 
D2B X,Y 2.888 14.438 0.0 
D2B X,Z -2.896 -14.481 0.0 
D2B Y,z 459 2.293 0.022 
D4 X,Y 1.783 11.959 0.0 

D4 X,Z -1.894 ~12.706 0.0 

D4 Y,Z -.695 -4.659 0.0 
DSA X,Y . 206 1.293 0.196 
DSA X,Z -.628 -3.948 0.0 
DSA Y,2 . 302 1.900 0.057 
DSB X,Y -.911 -5.904 0.0 
DSB Maz -.903 -5. 850 0.0 


DSB Y,Z -.203 -1.443 0.149 


7) 
w?| 


4. To get a significance of greater than .05 in this test 
requires a value between +/- 1.96, so only 3 of the 18 pairs 
of correlation coefficients are statistically equal, giving 
very strong evidence that the correlation coefficients of 
the residuals as recorded by the two sensors in a crossover 
Pair are, in general, not equal. 

The practical significance of the preceding tests is 
that there appears to be much local variation in the range 
and that one single model for simulating random replications 
of underwater track 1S not apparent. Therefore, each case 
must be simulated and studied separately, and the study of 
the error estimates distributions done on a case by case 


basis. 


IV. SIMULATION MODEL 


The Simulation model itself is a logical extension of 
the residual generation program. In fact, the method used 
to simulate a data track was to simulate residuals and = add 
them to the fitted straight line obtained for that track, 
rather than to start from scratch and simulate a totally new 
track at each iteration. This method was very quick and 
yielded very good simulated track segments. Figure 4 on the 
following page is graphical comparison of a typical original 
track segment overlaid by a track segment simulated by the 
method stated above. The comparison demonstrates the 
realistic quality of the simulated track. 

The regression method used to compute residuals of a 
track segment also produced the best straight line in three 
dimensional space to approximate the target path through the 
water. Given the straight line, the residuals themselves, 
and the assumption of normality of the residuals, 1t 1s 
possible to simulate a set of residuals from normal (Om ) 
deviates and add them to the straight line to obtain a 
Simulated track. 

In simulating residuals, the object is to compute an 
Nx 3 matrix, using normal (0,1) deviates, whose vectors of 
X, YY, and Z components are normally distributed with mean O 


and whose covariance matrix is equal to the covariance 
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matrix of residuals. That is 


Read x x’ 
where 
R= 3 x N matrix of simulated residuals 
XK = N x 3 matrix of N(0O,1) deviates 
T= 3 x 3 transformation matrix such that 


T x T’ = covariance matrix of residuals. 


UN 


Note that R must be transposed to get the desired N_ x 
residual format. 
It iS easy to show that any 3 x 3 matrix T will not 


alter the mean of QO. 


EGRI = ELT «x @ a 
= < Exe! 
Since the expectation of X -~ consisting of normal (0,1) 
deviates -- is identically O, 


SER] © Texte (0) ge O 
Giving the desired result. 
The condition that T x T’ is equal to the covariance 
matrix of the residuals is necessary because 


COVER] ECR x R’I/N (the mean is OQ) 


ELT x X¥’ «x X x» T’IZN 


T « ECX’ x XI/N x T’ 


(since T is a linear operator) 
— as Ga eae ees 
(because the covariance matrix of N(0O,1) random variables is 


the identity matrix) 
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The problem now becomes one of finding a 3 xX 3 matrix such 
that 

COV[IR] = T x T' 
We have the covariance matrix of the residuals. Let T be an 


upper triangular matrix. Then 


es Ti2 Tix Ls 0) 0) COV, COV, 5 COV, = 
.) T559 To< x Th9 T59 0) = COV... COV..5 COV... 
on. Sees "3 [es "Ss| gies eee aes 
Noting that the covariance matrix 1S symmetric and 


performing the matrix multiplication yields 








To< = SORT ‘COV_.) 
geese 
235 = 
33 
a an 2 
T59 = SORT (COV,5 T= ) 
COV,= 
T = 
13 t 
Ss 
aes Vee 
T = 
12 - 
22 
~ _ Zee 2 
bom = SQRT (COV, , Tix Th5 ) 


Since the matrix T is easily computed and X can be generated 
from a random number generating package, the residuals, R, 
are quickly and economically computed for as many simulated 


tracks as desired. 
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Adding the residuals thus formed to the best straight 
line target path of the original data yields a simulated 
track. 

The commented FORTRAN program written to perform the 
Simulation is included as Appendix A. It is basically a 
driver program that generates a user specified number of 
simulated tracks and interfaces with KEYMAIN, which 
estimates displacement and angular rotation values. Because 
some of the user friendly attributes of KEYMAIN interfere 
with the speedy generation of the simulation’s required 
output parameters, the package has been altered somewhat to 
increase speed. The output of the driver program is two 
data files. One data file, ROTATE.DAT, is an WN x 4 matrix 
that gives, in the first column, the maximum angle of 
rotation of the sensor, and, in the last 3S columns, provides 
the ordered Euler angle components that form the single 
maximum rotation angle. The second file, DISPLACE.DAT, is 
also an N x 4 matrix. The first column is the magnitude of 


displacement of the array, while the last S$ columns give the 


X, Y, and Z axis components of displacement. 


V. INTERFACE WITH EXISTING PROGRAMS 


The software developed to compute residuals and generate 
Simulated track segments are the original work of the 
author, aided by such canned routines as eigensystem 
analysis, random number generators and vector arithmetic. 
These canned subroutines came from the IMSL Library of 
Mathematic and Statistical functions developed for the IBM 
PC computers CRef.4], and are compiled in machine’ language 
libraries not reproducible here. The author ’s FORTRAN code 
is listed as Appendix A. 

In order to produce the estimates of sensor displacement 
and rotation, however, it was necessary to interface with a 
large stand alone FORTRAN package, KEYMAIN. This'7 program 
takes as input crossover region data sets (real or 
Simulated) and produces the estimates of displacement = and 
rotation of the second sensor of the crossover pair, based 
on the assumed accurate position of the first. The 
Orientation correction output by KEYMAIN is actually a three 
valued vector of ordered angular rotations in the XY, XZ and 
YZ planes. Similarly, the location correction is a three 
Valued vector of displacements in the X, Y and Zz directions. 
To reduce complexity, the ordered Euler angles were reduced 
to a single maximum angle of rotation about an appropriately 


tilted axis, and the three components of displacement were 
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reduced to a single quantity, magnitude of displacement. 
Thus the six dimensional quality of position corrections was 
reduced to two. 

KEYMAIN was designed as a stand alone product and as 
such 185 very user friendly. It was also designed to take 
not just one, but several, crossover data sets and produce 
displacement and rotation estimates for several sensors at 
once. Because its use in the simulation was to process a 
single simulated crossover data set, and because it was 
being called as a subroutine rather than used as a stand 
alone package, Significant changes were required in 
the program to obtain fast simulated results uninterrupted 
by the now unnecessary user friendliness. Not only did this 
require the modification of the executive driver routine, 


but also modification of several of the called subroutines. 
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VI. RESULTS 


It is imperative that the simulated track segments 
“look" and "act" like the original track segment they are 
Simulating. That the simulated track segments act like the 
Original is guaranteed by the equations: the residuals will 
have mean zero by definition and their covariance matrices 
were computed to be the same as the original's. For visual 
verification, Appendix D is included. Appendix D contains 
plots of each original track segment overlaid by one of the 
tracks simulated from it. A good fit 1S apparent in all 
cases. 

The simulation model can produce the data needed to 
produce a scatterplot of magnitude of displacement (in feet) 
versus maximum angle of rotation (in radians) like the 
schematic shown in Figure 5S. This represents an imaginary 
Situation where a target vehicle was driven through a 
Single crossover region on arange 700 times over the exact 
same track, and the results were fed into the program to 
produce the displacement and rotation values. It acts as a 
data base, or sampling distribution, for an array whose 
position has not changed, and the graph depicts the natural 
variability of the data. The contours represent some 
theoretical confidence levels for that particular sensor. 


For example, take the outermost contour, labeled .90. In a 
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future tracking exercise, crossover data can be fed into the 
program that produces displacement and rotation estimates. 
That will produce a point on the graph. If that point lies 
outside the contour, we have good evidence that the sensor 
has moved. Statistically speaking, one would have a 10 
percent chance of rejecting a true null hypothesis, givena 
null hypothesis of no sensor movement. Conversely, if the 
point plotted closer to the middle of the graph, there would 
be insufficient evidence to support the hypothesis of sensor 
movement, and the apparent movement would be attributed to 
the inherent variability in the data. 

It is unreasonable to expect that a target vehicle could 
be driven along the same track 700 times, nor would the 
range operators be likely to attempt it. The simulation 
program, however, needs only one straight line segment of 
track through a crossover region, and can then generate as 
many track data sets as needed to produce the graph. 

Figures 6 through 11 are plots produced from the 
Simulation model using the data sets from Figure 3 (p. 18). 
Figure 9 is a "well-behaved" plot that could conceivably, 
with many more runs, yield the type of graph displayed in 
Figure 5. In fact, the points appear so evenly distributed 
that one could conceivably assume independence between the 
rotation and displacement values. 

Although computationally attractive, independence is 


thought to be a poor assumption, because if a sensor fas 
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displaced a great distance, it has also had ample 
opportunity to rotate. This observation is speculative, but 
seems to be borne out in the remaining figures where a 
strong positive correlation appears to exist between the 
displacement and rotation values. The fishhook appearance 
on these graphs is aresult of making the correction 
estimates positive regardless of the direction of 
displacement’ or rotation. The graphs taken as a whole serve 
to illustrate the variable nature of the track data. 

These three graphs, drawn from the simulation, vividly 
illustrate the need for further study. Two points 
concerning the data on which these graphs were based are 
perhaps of some importance and may begin to explain the 
strange nature of the graphs produced. First, the data is 
fairly old, taken in September, 1982. The range operators 
from NUWES state that their capabilities have improved in 
the interim 3 years and that newer data could prove 
Significantly more accurate. Second, the data was taken on 
a single day, from a single range, using only 3 of the 24 
sensors on the range. 

The great variety of graphs produced points to the 
necessity to do much more research and to the utility of a 
Simulation model to accomplish it. In addition to the 
aforementioned use of newer data, the following 
considerations warrant study to ascertain their possible 


effects: 


ol 


(a) 


(b) 


(c) 


(d) 


(e) 


(f) 


Day to day variations on the range -- It is 
not clear at this time whether the character of 
the water column in which the target operates 
is invariant over time. Differences in salinity 
and/or temperature could conceivably evolve over 
time, affecting the accuracy of the track data. 


Seasonal variations on the range -- While 7 ie 
is not clear whether changes occur on the range on 
a daily basis, it certainly seems reasonable to 
expect a variation from season to season. Our 
data, taken from a single range on a single day, 
was insufficient to explore this. 


Location on the range-- Current practice on a 
tracking day is to take one sample of the water 
column on the range, checking for salinity, 
temperature, and other factors that affect the 
sound velocity profile, and assume that the 
results hold true for the duration of the exercise 
for the entirety of the range. The possibility of a 
daily change was discussed earlier. Here we mention 
the possibility of different water characteristics 
from one end of the range to the other. 


Geometry of tracking runs -—- It is possible 
that varying the geometry of the tracking run 
could result in changes to the quality of the 
data recorded. It can be seen from figure 3 (on 


page 18) that all of the data used (for the 
simulations was from tracks that run predominantly 
downrange with relatively little crossrange change 
and virtually no depth change. (Although depth 
cannot be seen from figure 3, it was examined for all 
the tracks.) 


Depth of target -- Depth appears to be the least 
reliable of the three dimensions recorded during 
tracking runs. Our data is from targets that 
operated in a narrow depth band. Deeper or 


Shallower targets could significantly affect the 
data. 


Inhomogeneities in the water -- It is quite 
conceivable that undetected inhomogeneities in the 
water could significantly affect the quality of 
the data. Tt would appear that currents and 
turbulence could affect the passage of sound 
through the water column, but quantifying these 
disturbances could be a major practical problem. 


a) 
tJ 


Precisely Characterizing the variability of the 
correction estimates is elusive and the preceding 
considerations invite much future research before that goal 
is reached. The existence of this simulation model 


brightens the outlook for ultimate success. 


CA 
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VII. AREAS FOR FUTURE WORK 


Problems encountered in the development of the 
Simulation model and intuition gained as a result of working 
on it gave birth to several areas for potential future 


study. Some of these are listed below. 


A. RESIDUALS FOR CURVED CROSSOVER DATA TRACKS 

It was convenient in this simulation model to use 
Straight line segments of track to get residuals. This 1s 
because the method of principle components works only for 
straight lines in N-space, rather than for curves. If a 
method could be developed to regress the data for any track 
onto its best path, regardless of curvature, one major 
obstacle in the use of the simulation would be removed. 
Whereas now we are limited in the type of data we can use, 


such a method would enable the use of all crossover data. 


B. WEIGHTED DATA BASED ON DISTANCE TO SENSOR ARRAY 

It was assumed in the model that recorded track data 
was uniformly accurate and reliable, regardless of the 
distance from target to sensor. That is, no distinction was 
made between the accuracy of the data when a target vehicle 
was "close" to a sensor and when the target was far away. 
Some sort of weighting scheme may prove beneficial and help 


smooth out the current data discrepancies. 
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C. TRISENSOR CROSSOVER DATA 

Although smaller in size, there exist areas on the 
range where a target may be tracked simultaneously by three 
sensors at once. These trisensor crossover regions could 
provide insight into the accuracy of the position of a 
sensor and may yield more accurate estimates of displacement 


and rotation. 


D. TRACK DATA SELECTION PROCEDURES 

It 15S mot now known whether a greater number of data 
points in a track segment yields better results in the 
Simulation and intuition is limited on this point. Research 
in this area could provide valuable guidelines for selection 


of data to use in the simulation in the future. 


E. DECISION RULES TO DETERMINE SENSOR MOVEMENT 

After accurate characterization of the variability of 
track data and correction estimates is made, the next 
logical and extremely useful step would be the construction 
of statistically sound decision rules to determine the 
following: 

(a) The discrepancy noted between the tracks in a 
crossover data set can be explained by the inherent 
variability of the data. 

(b) The discrepancy noted between the tracks in a 
crossover data set can be quantified by the sensor 
slippage model and can be computationally corrected. 

(c) The discrepancy noted between the tracks in a 
crossover data set cannot be explained by the model 


of sensor movement and some other explanation must be 
sought. 


Clearly, the last option is the least desirable, but if 


that situation is present, it needs to be noted. 


F. COMPUTATIONAL RANGE RESURVEY 

If the method of providing correction estimates can be 
proven to be reliable and accurate, it provides a potential 
method of range resurvey that has several advantages over 
the current method. First, the range would not need to be 
shut down to resurvey. In fact, range use would be 
mandatory to keep up to date sensor array positions. 
Second, it would be less expensive than the current method, 
replacing the equipment and manpower intensive current 
process with relatively inexpensive computer assets. Third, 
it would take less time. Resurvey of a single sensor array 
on the range can take up to a day; generating corrections 
from track data takes seconds. Fourth, since the current 
method uses a craft on the surface equipped with a pinger, 
all the pings must travel through the first 150-200 feet of 
the water column, where the sound velocity profile is quite 
variable and most difficult to determine, to get to the 
sensor array. In contrast, the underwater target vehicles 
tracked by the arrays are typically in 400-600 feet of water 
where the sound velocity profile is much smoother and easier 
to predict. Thus, one source of variability in determining 
sensor position is reduced. 

Tt is not envisioned that this computer method could 
ever replace the current survey process, but rather augment 
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it. Some way to determine the position of one sensor is 
necessary before KEYMAIN can even begin to function. 
However, if this computer process could augment current 
resurvey efforts, Or reduce the frequency with which 
resurveys must be made, significant savings in time and 


money could result. 


a 


APPENDIX A 


FORTRAN LISTING FOR PROGRAM SIMDAT 


PROGRAM SIMDAT2 


--This program simulates 3-D track data based on a 
-specific real track specified by the user. 
-User inputs are: 


one 1. track segment data file to be simulated 

aie 2. number of the left and right sensor in the 

os crossover pair 

ane 3S. mumber of simulated tracks desired 

one 4. request for sample simulated track (YES or NQ) 
as Je random number generator seed 


-The program output iss: 


1. 


file of residuals from the original track 
(RESIDUAL. DAT) 
file of a simulated crossover track, if 
requested (SIMTRACK. DAT) 
file of displacement values (in feet) for the 
right sensor of each simulated track in 4 
columns 

Gol 1i magnitude of displacement 

Col 2-4 X,Y,Z components of displacement 
file of rotation values (1n radians) for the 
right sensor of each simulated track in 4 
columns 

Col 1 maximum angle of rotation 

Col 2-4 ordered Euler angles of rotation 


-VARTABIE DECLARATION 


+ 


INTEGER*#4 N, I, J, K, IER, BIG1(3), BIG2Z(3), SIMS, 
POINT(130), TRACKS, IDL, IDR, TRKOUT 


CHARACTER DSNAME#13 


REAL *4 


NORM (260) 


REAL*#8 TRACK(130,6), TR1(130,3), TR2(130,3), MBARI (3), 
MBAR2 (3), MBARM1(130,3), MBARM2(130,3), COV1(4,3), 
COV2(3,3), DIF1, P1i(3,3), MUT1(130,3), SIM1(3,130), 
DD1¢(3), DD2(3), PA(3,3), PB(3,3), WORK(130), Di, D2, 
DIF2, SUM1, SUM2, Al. AZ. Ziteetaoyee zo (eptooe 
P2(3,35), ZT1(130,3), ZT2(130,3), TBAR, Z1iBAR, Z2EAR, 
TRKSUM(130), TKSUM2, CT1(3,130), CT2(3,130), SEED, 
MUT2(130,35), RESID1(130,3), RESID2(130,3), SU2(3,4), 


t++t+tettet 


38 


+ COV1R(3,3), COV2R(3,3), SQ1(3,3), SIMTRK(200,4), 
+ SIM2(3,130), ROTATE(1000,4), DISP(1000,4), DAIA(2,4) 


BB 
C...-Begin the user input section 
C 
WRITE (#,*) ‘Enter FNAME.FT of crossover data set’ 
WRITE (#.%) ° on disk : 
READ (*,°¢(A)’) DSNAME 
WRITE (*, *) 
WRITE (*,%*) ‘What is the NUMBER of the left sensor in’ 
WRITE (*.*) ° the crossover pair ? 
READ (*,*) IDL 
WRITE (#, *) 
WRITE (*#,*) ‘What is the NUMBER of the right sensor ? 
READ(*,*) IDR 
WRITE C#,*) °° 
WRITE (*,*) “How many simulated tracks do you desire ' 
WRITE (*#.%) ‘(NOTE = max 1000) ‘ 
READ (*,*) TRACKS 
WRITE (#,%) ° *° 
WRITE (#,*) “Do you want a sample simulated track ?' 
WRITE (*#,*) ‘Enter 1 for YES, O (zero) for NO 
READ (*,*) TRKQUT 
WRITE (#,*) ° * 
WRITE (*#,*) ‘Seed for the random number generator - 
WRITE (*#,%) “NOTE =: Seed must include a decimal 
READ (*,*) SEED 
WRITE (*, *) 
& 
C...-Read data from file 
S 


OPEN (1,FILE=DSNAME,STATUS=‘OLD‘) 
N = 0 
10 N= N + 1 
READ (1,*,END=30, ERR=30) POINT(N) , (TRACK (N, I) ,I1=1,5) 
B 
C...-Separate into "left" and "right" sensor tracks 
& 
DO 20 I = 1,3 


TRICN,I) = TRACK (N, I) 
TR2(N,1I) = TRACK (N, I+3) 
20 CONTINUE 
GOTO 190 
30 CLOSE (UNIT = 1) 
N= N- 1 
e 
C...-Compute the covariance matrix for each track. 
C...-First step, get column averages (with first column 
C... average, TBAR, computed for later use) 
C 
TBAR = O. 


po SO I = 1,3 


09 


MBARI(I) = O. 
MBAR2(I) = 0. 
DO 40 J = 1,N 


MBAR 1 (1) MBAR1 (I) + TR1I(J,1) 


MBAR2(T) MBAR2(I) + TR2¢(J,1) 
IF (I .EQ@. 1) TBAR = TBAR + DBLE(POINI (J)) 
40 CONTINUE 
MBAR1 (1) MBAR1 (I) / DBLE(N) 


MBAR2Z(I) = MBAR2(1I) / DBLE(N) 
70 CONTINUE 
TBAR = TBAR 


~N 


DBLE (N) 

Cc 

C..e-Do the matrix multiplication =: A(t)xA , subtracting oft 
C..e-the mean from each column entry to form the covariance 


C.eeematrix. Also make the matrix of (points —- means) for 
C...-iater use. 
cS 


DO 80 I = 1, 
DO 70 J = 1, 
COV1(I,J) = 
COV2(I,J) = 
DO 60 K = 1 


0. 
OQ. 
JN 
COV1(I,d) = COV1(1,3)+(TR1(K,I) -MBAR1 (I) ) 


* *(TR1(K,J) -MBAR1 (J) ) 
COV2(1,J) = COV2(1,3) +(TR2(K,1) -MBARZ(I)) 
* *(TR2(K,J) -MBAR2 (J) ) 
50 CONT INUE 
COV1(I,J) = COVI(I,J) / DBLE(N-1) 
COV2(1,J) = COV2(I,J) / DBLE(N-1) 
70 CONTINUE 


80 CONTINUE 
bO 100 I = 1,N 
DO 90 J = 1,3 
MBARM1 (I,J) 
MBARM2 (I,J) 


TRI¢(I,J) —- MBARI (J) 
TR2¢1,3) - MBARZ (J) 


H ll 


90 CONT INUE 

100 CONTINUE 
© 
C..-Form the matrix FP of ordered principle components for 
C...-@ach track. The columns of FP are the eigenvectors 
C...-associated with the eigen-values of the covariance 
C..ematrix for each track arranged in order of descending 
C...-e€igenvalues. (ie. the eigenvector associated with 
C...the largest eigenvalue is the first column) 
C 


C..-Call routine to compute eigenvalues/vectors for each 
Coes thace 


Cc 
CALL EIGRS(COV1,3,11,DD1,PA,3,WORK, IER) 
CALL EIGRS(COV2Z2,3,11,DD2,PB,3,WORK, IER) 
Cc 
C...Get eigenvectors in eigenvalue order, largest to 
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C...smallest, 


by column 


e 
C...-DDi/2 = eigenvalue vector 
C..e-PA/PB = eigenvector matrices 
Cc 
BIGi(1) = 1 
BIG2Z(1) = 1 
BIG1I(3) = 3 
BIG2(3) = 3 
IF (DD1(BIGi(1)) .LT. DD1(2)) BIGiI(1) = 2 
IF (DD2(BIG2(1)) .LT. DD2¢2)) BIG2(1i) = 2 
IF (DD1I(BIGiI(1)) .~.LT. DD1I(3)) BIGiI(1) = 3 
IF (DD2(BIG2(1)) .LT. DD2¢(3)) BIG2Z2(1) = 3 
IF (DD1¢(BIG1¢(3)) .GT. DDi¢1i)) BIGi(3S) = 1 
IF (DD2(BIG2¢(3)) .GT. DD2¢(1)) BIG2(3) = 1 
IF (DD1I¢(BIG1(3)) .GT. DDi(2)) BIGiI(3) = 2 
IF (DD2(BIG2(3)) .GT. DD2¢(2)) BIG2(3) = 2 
IF ((BIGi(1) + BIGi(3)) .E@. 3) THEN 
BIGi(2) = 3 
ESE 
IF ((BIGi¢1) + BIGi(3)) .EQ. 4) THEN 
BIGi¢(2) = 2 
EESE 
BIGi1(2) = 1 
END IF 
END IF 
IF ((BIG2¢(1) + BIG2¢(3)) .EQ. 3) THEN 
BIG2(2) = 3 
SiS] 3 
IF ((BIG2(1) + BIG2(3)) .~.EQ@. 4) THEN 
BIGZ(2) = 2 
ELSE 
BIG2(2) = i 
END IF 
END IF 
DO 130 I = 1,3 
DOD 120 J = 1,3 
P1(I,J) = PACI,BIGI(J)) 
P2(I,J) = PB(I,BIG2(J)) 
120 CONT INUE 
130 CONTINUE 
& 


C...-Compute the matrix ZT for each track 
C...ZT represents the projection of the track data onto 
C..e-the principle components 


ete ZU (TR — MBAR) 
C...-where TR — MBAR 
C...for each row 

C 


xP MBARM x P 
track data minus the column average 


C...Call routine to multiply matrices AxB 


C 


CALL VMULFF (MBARM1,P1,N,3,3,130,3,ZT1,130, IER) 
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CALL VMULF 
B. 
C...-Since there 
C...perform a si 
C...the first pr 
Cc 
TKSUM2 = O 
ZiBAR = Oo. 
Z2BAR = QO. 
DO 140 I = 
ZiBAR = 
Z2BAR = 
TRKSUM ¢ 
TKSUM2 
140 CONTINUE 


F (MBARM2,P2,N,3,3,130,3,Z212,130, IER) 


are some points missing from the data set, 
mple least squares linear regression onto 
inciple component 


1,N 
ZiBAR + ZT1(I,1) 
Z2BAR + ZT2(1,1) 
I) = DBLE(POINT(I)) -— TBAR 
= TKSUM2 + TRKSUM(I) #*2 


ZiBAR = Z1iBAR / DBLE(N) 


ZZBAR Zz 


u) 
e 
a 4 
be 
i wou Wl 


150 CONTINUE 
SUM1 
SUM2 
ZiBAR 
ZZ2BAR 
DO i160 [I = 
Zeist i 
Zee GL sd 
160 CONTINUE 
C 
C...Get CT matri 
C...0f the data 
C...-principle co 
& 


D 
ros 
lo uw ul Wt 


DO 170 I = 
Zi¢i,I) 
Z2¢1,1) 
170 CONTINUE 
DO i180 I = 
Z1¢<(2,1) 
Z1 (3.1) 
Levee td 
22S a1.) 
180 CONTINUE 
Cc 
C...-Call routine 
Ee 
CALL VMULF 


BAR / DBLE(N) 


1,N 

(ZT1¢I,1) -— ZiIBAR) * TRKSUM(I) 
(ZT2(I,1) -— Z2BAR) * TRKSUM(I) 
SUM1i1 + DIFi 
SUM2 + DIF2 


/ TKSUM2 
/ TKSUM2 
—- Di * TBAR 
—- D2 * TBAR 
1,N 
) = AL + D1 * DBLE(POINT(I)) 
) = A2 + D2 * DBLE(POINT(I)) 


x which represents the orthogonal projection 
onto the straight line of the first 
mponent 


1 | 


3 N 
ZT1i¢1,1) 
ZT2¢1,1) 


COoOCOo2 


ooo 0 


| | | 


to multiply matrices AxB 
an Ga? |e Ae UE: San | FS He BG fo RS YS) SS) 


62 


CALL VMULFF (P2,22,3,3,N,3,3,CT2,3, IER) 
Cc 
C..-Move "line" of data back into original coordinate system 
Cc 
DO 200 I = 1,N 
DO 190 J = 1,3 
MUT1(I,J) = CT1(J,I) + MBARI (J) 
MUT2¢I1,J) = CT2(J,1) + MBAR2(J) 


Cc 
C...-Compute residuals for each track 
C 
RESID1I(I,J) = TRI(I,J) — MUT1I(I,J) 
RESID2(I,J) = TR2(I,J) -— MUT2(I,J) 
190 CONTINUE 
200 CONTINUE 
Cc 


C..e-Write the set of residuals out to the file RESIDUAL. DAT 
C 
OPEN (2, FILE = ‘RESIDUAL.DAT’, STATUS = ‘NEW’) 
DO 210 I = 1,N 
WRITE (2,330) (RESID1I (I,J) ,J=1,3) , (RESID2¢(1,J) ,J=1,35) 
210 CONTINUE 
C 
C...-Compute the covariance matrix of the residuals 
C...(Note =: column averages are identically zero) 
C 
C...Call routine to multiply matrices A(t)xB, then divide 
€..-by N-1 
C 
CALL VMULFM(RESID1 ,RESID1 ,N,3,3,130,130,COV1IR,3, 1ER) 
SAVE MULPM(RESID2,RESIDZ,N,S,5, 150, 150, COV2R3s, IER) 
DO 230 I = 1,3 
DO 220 3 = 1,5 
COVIR(I,J) 
COV2R (I,J) 


COVIR(I,J) / DBLE(N-1) 
COV2R(I,J) / DBLE(N-1) 


220 CONTINUE 
230 CONTINUE 
C 
(C...Get "Square root" of covariance matrix of residuals 
. 
SQ1(23,3) = DSORT(COVIR(S,S)) 
SQ2(¢3,3) = DSORT (COVER (S,3)) 
SQ1(2,3) = COVIR(2,3) / SQ1(3,35) 
SQ2(2,3) = COV2R(2,3) 7/7 SQ2(5,5) 
SQ1(2,2) = DSQRT(COVIR(2,2) — SQ1 (2,3) **¥2) 
SQ2(2,2) = DSORT(COV2R(2,2) - SQ2(2,35) **2) 
SQ1¢1,3) = COV1IR(1,3) 7/7 5SQ1¢(3,3) 
SQ2(1,3) = COV2R(1,3) / SQ2(3,3) 
SQ1(1,2) = (cCOViIR(1,2) — S$Q1(1,3)*S01(2,3)) / SQ1 (2,2) 
SQ2(1,2) = (COV2R(1,2) -— SQ2(1,3)*#5Q2(2,35)) / SQ2(2,2) 
SOT (1, 1) DSORT (COV1IR(1,1)—-(SQ1 (1,2) ¥¥2+SQ1 (1,5) **2) ) 


SQ2¢(1,1) DSORT (COV2R (1,1) -(SQ2¢61,2) **¥2+5Q2(1,5) **2)) 
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SQ1(2,1) ==0. 
SQ2(2,1) = 0. 
SQ1(3,2) = 0. 
SQ2(3,2) = Oo. 
SQ1(3,1) = 0. 
SQ2(3,1) = 0. 


& 
C...Compute sets of residuals and get rotation/displacement 
C...values 
S 
DO 410 SIMS = 1,TRACKS 
& 
C...Compute set of simulated residuals from normal (0,1) 
C...deviates 
C 
DG 260 I = 1,3 
& 
C...Call routine to generate Normal (0,1) deviates 
S 
CALL GGNPM(SEED,2*N,NORM) 
DO 2530 J = 1,N 
RESID1¢(J,1) 
RESID2¢(J,1) 
250 CONTINUE 
260 CONTINUE 


NORM (J) 
NORM (N+J) 


Bs 
C..2e-Call routine to multiply matrices AxB(t) 


C 
CALL VMULFP (SQ1,RESID1,3,3,N,3,130,SIM1,35, IER) 
CALL VMULFP (SQ2,RESID2,3,3,N,3,130,51IM2,3, I1ER) 

& 

C...Put together Nx6 matrix of simulated tracks for both 


C...arrays By adding straight line in original coordinate 
C...Ssystem to residuals 
fc 
DO 280 I = 1,N 
DO 270 J = 1,3 
SIMTRK(I,J) = SIM1(J,1I) + MUT1¢(I,J9) 
SIMTRK(I,J+3) = SIM2¢(J,1) + MUT2(1I,J) 
279 CONT INUE 
280 CONTINUE 
C 
C...Write the first simulated track out to the file 
C..-SIMTRACK. DAT if a sample simulated track was requested. 
eS 
IF((SIMS .£EQ. 1) .AND. (TRKOQUT .GT. 0)) THEN 
OPEN (3, FILE = ‘“SIMTRACK.DAT’, STATUS = ‘NEW’) 
DQ 285 I = 1,N 
WRITE (3,340) (SIMTRK(I,J),J3 = 1,6) 
293 CONTINUE 
END IF 
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C...Feed the simulated track into KEYMAIN to get rotation 
C...and displacement numbers 
Bi 
CALL KEYSUB(SIMTRK,N,DATA, IDL, IDR) 
C 
C...-Make 2 matrices - one for displacement data and one for 
C...the rotation data 
C 
DO 290 I = 1,4 
DISP(SIMS,1I) = DATA‘1,1) 
ROTATE(SIMS,I) = DATA(2,1) 
290 CONTINUE 
WRITE (#,320) SIMS 


S 
C...-60 back and do it again 
S 
410 CONTINUE 
& 


C..-After TRACKS simulated tracks, write the displacement 
C...sets and the rotation sets out to a file 
oS 
OPEN (4,FILE "DISPLACE.DAT’,STATUS = ‘NEW’) 
OPEN(S,FILE "ROTATE.DAT’,STATUS = ‘NEW’) 
DO 300 I = 1,TRACKS 
WRITE (4,310) (DISP(I,J),Jj = 
WRITE (S,310) (ROTATE(I,J) J 
300 CONTINUE 


1,4) 
= 1,4) 


& 

C..-Close out the files 

& 
CLOSE(CUNIT = 2) 
CLOSE(UNIT = 3) 
CLOSE (UNIT = 4) 
CLOSE(UNIT = 3) 
STOP 


310 FORMAT (2X,4F17.8) 
320 FORMAT(2X,’Through KEYMAIN ‘,146,’ time(s) so far‘) 
330 FORMAT(1X,6F12.7) 
340 FORMAT(1X,46F12.5) 
END 
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APPENDIX E 


FORTRAN SUBROUTINES CONECT AND REDUCE 


The program KEYMAIN requires that any array for which 
correction estimates are desired must be "connected" to the 
first input sensor array by one, or a series of, crossover 
data sets. For example, if an input crossover data set uses 
sensor arrays 3 and 6, while another uses arrays 6 and 10, 
all three of the arrays (3S, 6 and 10) are connected. Me ge 
however, a first input data set uses arrays 7 and 7%, a 
second input data set uses 12 and 10, while a third input 
data set uses 9 and 13, the arrays 7, %? and 13 are 
connected, 12 and 10 are connected, but all five arrays do 
not form a single connected set. In this case, if 7 was the 
first array input to the program, correction estimates could 
not be made for arrays 10 and 12. The subroutine CONECI 
checks to see that connectedness exists in the input data 
before KEYMAIN is allowed to continue. 

KEYMAIN allows 3 options if CONECT discovers that the 
arrays of the input data sets are not connected. One option 
is to quit, in which case the program terminates. Another 
option 1s to add more data so that all the arrays are 
connected. In the second example above, for instance, if a 
crossover data set using arrays 9 and 12 was input into the 


program, all 3 arrays would then be connected and KEYMAIN 
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would continue. A third option is to continue with KEYMAIN, 
but use only the first connected set, that is, all the 
arrays that are connected to the left array of the first 
input crossover data set. This option presents special 
problems because it requires a reduction of the data 
structures that have been built as each data set was input. 
The subroutine REDUCE does this data structure reduction and 
is called from KEYMAIN only when this third option is 
selected. 

KEYMAIN passes to CONECT two pieces of information. 
The first is the variable Ri oan represents the number of 
data sets that were input into KEYMAIN. The second piece of 
information is a2x R1 matrix, IND2, that contains the 
number of the left and right sensors of the Ri crossover 
data sets. Row 1 contains the left sensor numbers, row 2 
the right. CONECT performs its connectedness check by 
starting a variable length list that contains the = array 
numbers of those arrays that are connected. The list starts 


with only two entries, the left and right sensor arrays of 


the first crossover data set. These are connected, and the 
left sensor is the "root" to which all arrays”) should 
connect. Elements are added to the list by sequencing 


through the list from the beginning, and adding to the list 
any array numbers for IND2 that (1) are not yet on the list 
and (2) are connected by a crossover data set to an array 


that is already on the list. If, after sequencing through 
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the variable length list, there are arrays remaining in IND2 
that never were put on the list, those arrays were not a 
part of the first connected set. If all arrays in IND2 were 
ein on the list, all the arrays were connected. After 
the first list is exhausted,  CONECT repeats this procedure 
as many times as there are disjoint sets of arrays, starting 
each new list with the arrays of a data set not in any 
orevious set. CONECT informs the user of the individual 
sets of connected sets of arrays and raises a flag to alert 
the user if all sets were not connected. 

If all the arrays in the input crossover data sets were 
not connected, the user has three options to proceed, 
described earlier. If the chooses to continue using the 
first set of connected arrays, REDUCE is called to pare the 
data strauctures built up during the data input process. In 
particular, the Ri x 3 x 3 array CROSSA and the Ri x 6 
matrix mean need to be reduced to contain only those 
elements that correspond to the data sets in the first 
connected set. CROSSA, which has Ri "nages" of 3 x 3 
matrices of crossproduct deviations from the mean, needs to 
have those pages removed that correspond to every crossover 
data set in the original input not connected to the first 
array. MEAN, which has Ri rows of the colunm averages for 
each data set, needs to have those rows removed that 


correspond to data sets not connected to the first array. 
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The variable Ri itself will be reduced to reflect this 
smaller subset of connected array pairs. 

CONECT stores the information that indicates which data 
sets in IND2 are connected to the first array. This 
information is passed to REDUCE through KEYMAIN. REDUCE 
reforms the data structures to reflect the smaller number of 
data sets now being considered by KEYMAIN. 

CONECT also provides the matrix IND1, aK x Ril matrix 
that is used elsewhere in KEYMAIN. The R1 columns represent 
the crossover data sets input into KEYMAIN. The variable K 
represents the number of individual arrays in the Rl data 
sets, each row representing a separate array. For each 
column in INDi, the entries are all O except in the rows 
that represent the left and right sensor for that column's 
data set. If the row corresponds to the left array, the 
Value is l. If it corresponds to the right array, the value 
is 2. If REDUCE is called, some columns (representing 
crossover data sets) and rows (representing individual 
sensors) of IND1 need to be omitted. REDUCE reforms IND1 to 


its smaller size. 
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APPENDIX F 


FORTRAN LISTING FOR SUBROUTINE CONECT 


SUBROUTINE CONECT (OUT,R1,IND2,K,IND1,IA,TESTC,IND2R, 


a DATSET) 
CG 
CE EE IE IE EE IE FE HE SEE IEE SEE HE IEE EE EE EEE EE HEHE IE HEE EEE EE HE EEE EH HEHE EE HE 
C This subroutine checks for the connectedness of the © 
Ey input data sets. If the problem is connected then the 
e user 1s informed and the array pairs are printed on the 
= screen; if not connected, then the user is prompted to 
Ee select one of three options —- quit, add conecting data 
C sets, or run the program using the first connected set 
C that was input. Gygax - July 1985 
CHEE HEHEHE HEE HEHEHE HEE EERE EEE EEE EEK EHR HEH 
Cc 
Cc ...-Variable declarations. 
C 
INTEGER*#4 R1I,K,IND2(2,30) , IND1 (30,30) ,I1,J,1A(30) ,FIRSI 
INTEGER#4 LIST(30) ,BEGIN,HALT,DISCON,L,M,O,TESIC,OU!, 
INTEGER#4 DATSET (30) , COUNT , SAVE (2,30), IND2R(2, 30) 
C 
C ---Initialize the values of FIRST and COUNT: 
C 
FIRST = 0 
COUNT = 9 
Us 
C ---Make vector [A = list of all arrays (w/o repeats) 
GS. in IND2 and get the value for K = # of individual 
C arrays. 
& 
TA(1) = IND2(1,1) 
TAC2) = IND2(2,1) 
K = 3 
IF (Ri .EQ@. 1) GOTO 690 
DO SO I = 1,R1 
DO 40 J = 1,2 
M= kK - 1 
DO SO L = 1,M 
IF (IND2(J,1) .EQ@. IA(L)) GOTO 40 
30 CONT INUE 
TACK) = IND2(J,T1) 
K =k + 1 
40 CONT INUE 
30 CONT INUE 
60 K=kK - il 
C 
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WRETE< GUT, *) Ri 


WRITE (OUT,*) 'K 


C «eeFor each column of IND! (columns correspond to 
C data sets) the entries are all zero except for 
C the row that corresponds to the left array ‘= 1) 
Ls and the right array (= 2). 
C 
DO 80 I = 1,Ri 
DO 7O 3g = 1,K. 
INDi(J,I) = O 
Peet tbeGt., 1) Mees sinc d)) INDICI,1) = 1 
De Orne (2 CER 1Akd)) INDICI,1I) = &2 
Tey CONT INUE 
BO CONTINUE 
Cc 
C ---Check to see if all the arrays are cennected. 
is 
MESTG = 1 
bro! (1) Sa=1AC1) 
DO 131 [ = 1,Rl1 
Peete tte. “May tly) INDe l,i) = -INDE(1, 1) 
Peo ween foe rieloiGlpesIND2Z(2et) = -INDZ(., 
131 CONTINUE 
BEGIN = 1 
HALT = 1 
Oe Gliese GtiN .LE. BALI) GOTO 1709 
NODE = LIST (BEGIN) 
REGIN = BEGIN + 1 
Catad ft = 15R1 
Pentecost Os PND? Clete AND. CIND2(2,f) .G7T. 0) 
# GOTO 150 
HALT = HALT + 1 
Pst ee — ID 2.(2. 1) 
DO 141 J = 1,R1 
Poo ei) ES KiekST (HALT) 6 INOS 1, J) =-TNDE (1D 
Pei 2). eels) CAL Tele LIND Gee, J) -— COMO | 
141 CONTINUE 
i3Q 8 =6CONTINUE 
wero Lose 1yRl 
eer Oem UNG Wee ei ate Ceo ky) AND. CINDS «1, ft) .o7. 
# E071UO so 
HALT = HALT + 1 
reomtnelt) = -—IND2(1, 1) 
Dome ted = bykl 
Pio te) ees —takol ety) ) END2¢1 ,J*=-TND? | 
fees (Pewee. -LiSl (Heise INS262,J)=-INOD 
Pol CONTINUE 
159 CONTINUE 
GOTO 140 
i CONTINUE 
5: 
C eePrink out the matched ovairs. 


‘Ri 
ty te 


. 
i 


% 


5 


170 


ZOO 


ans 


ye 


eet) 


2 
ae ows - 


240 


DISCON = O 

WRITE (OUT, 230) 

DO 200 I = 1,R1 

IF (IND2(1,I) .LT. 0) GOTO 190 
IF (IND2(1,1I) .EQ@. 0) GOTO 200 


TF (C€IND2¢1,1T) .GT. QO) .~. AND. (DISCON .EQ. 


FIRST = FIRST + 1 
DISCON = 1 

TESTC = 0 

BEGIN = 1 

HALT = 1 

LIST(1) = IND2(1,1) 

GOTO 200 

WRITE (OUT, 240) -IND2(1,1) ,-IND2(2,1) 


THEN 
COUNT = COUNT + 1 
IND2R(1,COUNT) = —IND2(1,1) 
IND2R(2,COUNT) = —IND2(2,1) 
DATSET (COUNT) = I 

END IF 

SAVE(1,1) = -IND2(1,1) 

SAVE (2,1) = -IND2(2,1) 

IND2(1,I1) = 0 

IND2(2,I1) = 0 

CONT INUE 


IF (DISCON .EQ. 1) GOTO 140 
DO 220 I = 1,R1 


IND2(1,1) = SAVE(1,1) 
IND2(2,1I) = SAVE(2,1) 
CONTINUE 
RETURN 


i) 


FORMAT CIX, THE FOLLOWING PATRS ARE CONNE@ he) 


FORMAT (1X,1415) 
END 


pe) 


GOTO 


ZOO 


IF ((CFIRST.E0.0) .GR. (CF IRST SES. 1) SAND. (DISEOMAay. 


APPENDIX G 


FORTRAN LISTING FOR SUBROUTINE REDUCE 


SUBROUTINE REDUCE (CROSSA,MEAN,R1,K,IND1,1A, 
IND2R,DATSET) 


HEHE HEHEHE HE HEH HEHEHE HEH HEHEHE EH EH HEHEHE HEHEHE EHH EHH EEE HH 


OO onIIan 


0 


O00 


O00 nm 


This is a specialized subroutine that is used when 
option three is invoked as a result of a failed con- 
nectedness test. The disconnected data sets must be 
removed from the variables CROSSA and MEAN, and other 
Program supporting variables must be adjusted. 

Gygax - July 1985 


HHH HEHEHE HEHE EEE EEE EEE EEE EERE EE HEE EEE HEE REE EHH HHH H 


10 
29 


-.. Variable declarations. 


INTEGER*4 R1,K,IND1 (30,30) ,IA(30) , IND2R(2,30) ,1,J,L, 
M,DATSET (30) 


REAL*8 CROSSA(30,3,3) ,MEAN(3S0,6) 
-.. Compute the new, reduced Ri: 


Da 10 Ila 1,30 
IF (IND2R(1,I1) .£Q. 0) GOTO 20 


CONTINUE 
Ri =1- 1 
»»- Make new, reduced vector IA = list of all arrays 


in IND2R w/o repeats. Also, compute a new EK. 


IA(1) 
IA(2) 
K = 3 
IF (Ri .EQ. 1) GOTO 40 
DO SO I = 1,R1 
DO 40 J = 1,2 
M=K-1 
DO 30 L = 1, 
IF (IND2R 
CONTINUE 
IA(K) = IND2R(J,I) 
K=K +1 
CONTINUE 
CONTINUE 
K=K - 1 


IND2R (1,1) 
IND2R (2,1) 


M 
(J,I) ~EQ. TA(L)) GOTO 40 
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Oooo n0g0ng 


a 


790 
80 


9O 


100 
110 
1290 


~2. Remake the reduced matrix INDI - for each column 
in IND1 (corressponding toa a data set) the 
entries are zero except for the entries corres-~ 
ponding te the left array (= 1) and the right 
array (= 2). 


DO 80 I = 1,R1 
DO 70 J = 1,K 
IND1(J,I) = 0 
IF (IND2R(1,1) -EQ. IA(J)) IND1I(J,I) 
IF (IND2R(2,1) .EQ@. IA(J)) IND1(J,1) 
CONTINUE 
CONTINUE 


i oil 
= 


».» Reduce the arrays CROSSA and MEAN to account 
for the removed data sets. 


DEMO t=) tert 
DO 90 J =.1,6 
MEAN(I,J) = MEAN(DATSET(I) ,J) 
CONT INUE 
rea] slo merep 5) 9S 
DO 100 L = 1,3 
CROSSA(1,J,L) = CROSSA(DATSET(I) ,J,L) 
CONT INUE 
CONT INUE 
CONTINUE 
RETURN 
END 
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